/*
 * Clase: MatrizDispersa
 *
 * Versi�n 1.0
 *
 * Copyrigth: Trabajo Pr�ctico Hector Morel - Gustavo Plan�s
 *
 */

package ecuacion;

/**
 * Clase Matriz Dispersa
 * @author Hector Morel, Gustavo Plan�s
 */

import java.util.Vector;

public class MatrizDispersa<E extends Numero> {
	
	/**
	 * Constructor
	 * @param Orden de la Matriz
	 */
	public MatrizDispersa(int x) {
		elementos = new Vector<E>(CAPACIDAD_POR_DEFECTO, CAPACIDAD_POR_DEFECTO); // Crea vectores
		filas = new int[CAPACIDAD_POR_DEFECTO];
		columnas = new int[CAPACIDAD_POR_DEFECTO];
		orden = x; // Orden de la Matriz
		n = 0; // Cantidad de elementos distintos de cero
	}

	/**
	 * Cargar Elemento
	 * @param Elemento, Fila, Columna en las que se desea cargar el Elemento
	 */
	public void cargarElemento(E elemento, int fila, int columna) throws Exception {
		
		int low = 0;
		int high = n - 1;
		int mid;
		int r = -1;
		boolean encontrado = false;

		// Verifica la inserci�n dentro de los l�mites de la matriz
		if (fila > orden || columna > orden || fila < 1 || columna < 1) 
			throw new Exception("Posicion fuera de rango");

		// Hace una b�squeda binaria para ver si existe un elemento cargado en esa posici�n
		while(low <= high)
		{
			mid = (low + high) / 2;

			if(filas[mid] < fila) {
				low = mid + 1;
				r = mid;
			} else if(filas[mid] > fila)
				high = mid - 1;
			else {
				if(columnas[mid] < columna) {
					low = mid + 1;
					r = mid;
				} else if(columnas[mid] > columna)
					high = mid - 1;
				else {
					elementos.setElementAt(elemento, mid);
					encontrado = true;
					break;
				}
			}
		}
		
		if (encontrado == false) { // Si no existe elemento en dicha posici�n, lo crea
			if (n + 1 == filas.length) // Si los vectores est�n llenos, los duplica
				duplicarVectores();

			if (r != n - 1) // Si no se inserta al final, se hace un hueco
				huecoVectores(r + 1);

			elementos.insertElementAt(elemento, r + 1); // Inserta en la posici�n indicada
			filas[r + 1] = fila;
			columnas[r + 1] = columna;
			n++;
		}
	}

	/**
	 * Suma una matriz a la matriz usada
	 * @param Matriz dispersa a sumar
	 */
	public void sumaMe(MatrizDispersa<E> mat) {

		int pos;
		if (mat.getOrden() == orden) { // Verifica que las matrices sean del mismo orden
			for(int i = 0; i < mat.n ; i++) {
				pos = binarySearch(mat.filas[i], mat.columnas[i]);
				if (pos > -1) // Suma
					elementos.setElementAt((E) elementos.elementAt(pos).sumar(mat.elementos.elementAt(i)), pos);
				else
					try { // Si no existe el elemento lo inserta
						cargarElemento(mat.elementos.elementAt(i), mat.filas[i], mat.columnas[i]);
					} catch(Exception e) {}
			}
		}
	}
	
	/**
	 * Suma dos Matrices Dispersas
	 * @param Matriz dispersa a sumar
	 * @return El resultado de la suma
	 */
	public MatrizDispersa suma(MatrizDispersa<E> mat) {
		
		int pos;
		MatrizDispersa<E> copia = this.copiar();
		if (mat.getOrden() == orden) { // Verifica que las matrices sean del mismo orden
			for(int i = 0; i < mat.n ; i++) {
				pos = copia.binarySearch(mat.filas[i], mat.columnas[i]);
				if (pos > -1) // Suma
					copia.elementos.setElementAt((E) copia.elementos.elementAt(pos).sumar(mat.elementos.elementAt(i)), pos);
				else
					try { // Si no existe el elemento lo inserta
						copia.cargarElemento(mat.elementos.elementAt(i), mat.filas[i], mat.columnas[i]);
					} catch(Exception e) {}
			}
		}
		return copia;
	}
	
	/**
	 * Halla la inversa de la matriz dispersa
	 * @return Inversa de la matriz
	 * @throws Exception
	 */
	public MatrizDispersa inversa() throws Exception {

		for(int i = 1; i <= orden; i++)
			existeFila(i); // Verifica que exista por lo menos un elemento distinto de cero en la fila; si esto no sucede, la matriz no es invertible
		MatrizDispersa<E> copia = this.copiar(); // crea una copia de la matriz original
		MatrizDispersa<E> identidad = new MatrizDispersa<E>(orden); // Crea una matriz identidad (Siempre Doble)
		E uno = elementos.elementAt(0);
		for(int i = 1; i <= orden; i++)
			identidad.cargarElemento((E) uno.devuelveUno(), i, i);
		E epivot; // Indica el elemento que se hatomado como pivot para dejar la mariz en forma canonica por filas
		int ppivot = 0; // Indica la posicion en la que se encuentra el elemento pivot
		int pfs; // Variable auxiliar que se utiliza para moverse en las filas de la copia en las que no se encuentra el pivot
		int pfsI;// Variable auxiliar que se utiliza para moverse en las filas de la Identidad en las que no se encuentra el pivot
		int pfaI = 0; // Variable auxiliar que se utiliza para moverse en la fila de la identidad en las que se encuentra el pivot
		int ppI = 0; // Variable auxiliar que indica donde empieza la fila en la que se ha tomado un pivot, en la identidad
		int pfa;   // Variable auxiliar que se utiliza para moverse en la fila de la copia en las que se encuentra el pivot
		int fa;  // Variable auxiliar que se utiliza para dar fin al ciclo que encuentra el primer elemento de la fila en la que se ha tomado pivot en la identidad
		E oper;
		int i;
			
		// PROCESO DE ESCALONAMIENTO
		for(i = 1; i < orden; i++) {
			copia.existeFila(i);
			int f  = i;
			while(copia.filas[ppivot] != i)
				ppivot++;// Selecciona un pivot
			epivot = (E) copia.elementos.elementAt(ppivot);
			pfs = ppivot;
			while(f == i)
				f = copia.filas[++pfs]; // Obtiene la fila siguiente a la del pivot
			while (identidad.filas[ppI] != i) // Hace lo mismo pero con la identidad
				ppI++;
			pfsI = pfaI;
			fa = i;;
			while(fa == i)
			fa = identidad.filas[++pfsI];

			while (pfs < copia.n) {
				pfa = ppivot;
				pfaI = ppI;
				f = copia.filas[pfs]; // Selecciona un operador para realizar las operaciones elementales sobre las filas
				oper = (E) copia.obtenerElemento(copia.filas[pfs], copia.columnas[ppivot]);
				if (!oper.esCero()) {
					while(pfs < copia.n && f == copia.filas[pfs] && copia.filas[pfa] == i) {
						if(copia.columnas[pfa] == copia.columnas[pfs]) {
							E mult = (E) copia.elementos.elementAt(pfs).multiplicar(epivot);
							copia.elementos.setElementAt((E) copia.elementos.elementAt(pfa).multiplicar(oper).restar(mult), pfs); // Si las columnas de los elementos que se encuentran en la fila que contiene al pivot y la que contiene al elemento que se desea dejar en cero son iguales, realiza la operacion elemental indicada
							if (copia.elementos.elementAt(pfs).esCero())
								copia.eliminarElem(pfs); // Si el elemento se ha cerado se elimina del vector de elementos
							else
								pfs++;
							pfa++;
						} else if (copia.columnas[pfa] > copia.columnas[pfs]) {
							copia.elementos.setElementAt((E) copia.elementos.elementAt(pfs).multiplicarInversa(epivot), pfs); // Si las columnas no coinciden, se asume que uno de los dos elementos es cero
							pfs++;
						} else{
							copia.cargarElemento((E) copia.elementos.elementAt(pfa).multiplicar(oper), f, copia.columnas[pfa]);
							pfa++;
							pfs++;
						}
					}
					// Si el ciclo no llego al final se fuerza a que esto ocurra
			 		if(f == copia.filas[pfs])
			 			while(f == copia.filas[pfs]) {
		 					copia.elementos.setElementAt((E) copia.elementos.elementAt(pfs).multiplicarInversa(epivot), pfs);
			 				pfs++;
						}
			 		if(copia.filas[pfa] == i)
			 			while(copia.filas[pfa] == i){
	 						copia.cargarElemento((E) copia.elementos.elementAt(pfa).multiplicar(oper), f, copia.columnas[pfa]);
			 				pfa++;
			 			}
			 		// Se hace el mismo procedimiento pero con la matriz identidad
					while(pfsI < identidad.n && f == identidad.filas[pfsI] && identidad.filas[pfaI] == i) {
						if(identidad.columnas[pfaI] == identidad.columnas[pfsI]) {
							E mult = (E) identidad.elementos.elementAt(pfsI).multiplicar(epivot); 
							identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfaI).multiplicar(oper).restar(mult), pfsI);
							if (identidad.elementos.elementAt(pfsI).esCero())
								identidad.eliminarElem(pfsI);
							else
								pfsI++;
							pfaI++;
						} else if (identidad.columnas[pfaI] > identidad.columnas[pfsI]) {
							identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfsI).multiplicarInversa(epivot), pfsI);
							pfsI++;
						} else {
							identidad.cargarElemento((E) identidad.elementos.elementAt(pfaI).multiplicar(oper), f, identidad.columnas[pfaI]);
							pfaI++;
							pfsI++;
						}
					}
					if(f == identidad.filas[pfsI])
						while(f == identidad.filas[pfsI]) {
							identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfsI).multiplicarInversa(epivot), pfsI);
							pfsI++;
						}
					if(identidad.filas[pfaI] == i)
						while(identidad.filas[pfaI] == i) {
							identidad.cargarElemento((E) identidad.elementos.elementAt(pfaI).multiplicar(oper), f, identidad.columnas[pfaI]);
							pfaI++;
						}
				} else
					while(pfs < copia.n && f == copia.filas[pfs])
						pfs++;
			}
		}
		
		// PROCESO DE RETROSUSTITUCION
		/* Al hacer la retrosustitucion la matriz ya debe estar en forma triangular y lo que se hace es repetir el proceso anterior pero en forma inversa,
		  es decir, se empieza por atras, y el ultimo elemento distinto de cero de cada fila sera el pivot, y se realizaran operaciones elementales sobre
		  las filas que se encuentren antes de fila en la que se encuentre el mismo
		 */
		pfa = copia.n - 1;
		ppivot = copia.n - 1;
		pfaI = identidad.n -  1;
		ppI = identidad.n  - 1;
		int f = 0;
		for(i = orden; i > 1; i--) {
			copia.existeFila(i);
			f  = i;
			ppivot = pfa;
			while(copia.filas[ppivot] != i) {
				ppivot--;
			}
			epivot = (E) copia.elementos.elementAt(ppivot);
			pfs = ppivot;
			ppI = pfaI;
			while (identidad.filas[ppI] != i) // hace lo mismo pero con la identidad
				ppI--;
			while(f == i) {
					copia.elementos.setElementAt((E) copia.elementos.elementAt(pfs).dividir(epivot), pfs);
					--pfs;
					if (pfs < 0)
						break;
					else
						f = copia.filas[pfs];
			}
			pfsI = pfaI;
			fa = i;
			while(fa == i) {
				identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfsI).dividir(epivot), pfsI);
				--pfsI;
				if(pfsI < 0)
					break;
				else
					fa = identidad.filas[pfsI];
			}
			if (pfs < 0 || pfsI <0)
				break;
			oper = elementos.elementAt(0);
			epivot = (E) oper.devuelveUno();
			while (pfs >= 0) {
				pfa = ppivot;
				pfaI = ppI;
				f = copia.filas[pfs];
				oper = (E) copia.obtenerElemento(copia.filas[pfs], copia.columnas[ppivot]);
				if (!oper.esCero()) {
					while(f == copia.filas[pfs] && copia.filas[pfa] == i) {
						if(copia.columnas[pfa] == copia.columnas[pfs]) {
							E mult = (E) copia.elementos.elementAt(pfs).multiplicar(epivot);
							copia.elementos.setElementAt((E) copia.elementos.elementAt(pfa).multiplicar(oper).restar(mult), pfs);
							if (copia.elementos.elementAt(pfs).esCero()) {
								copia.eliminarElem(pfs);
								pfa--;
								ppivot--;
							}
							pfs--;
							pfa--;
							if (pfs < 0)
								break;

				        } else {
							copia.elementos.setElementAt((E) copia.elementos.elementAt(pfs).multiplicarInversa(epivot), pfs);
							pfs--;
							if (pfs < 0)
								break;
						}


					}
				   if(pfs >= 0)
			 			while(f == copia.filas[pfs]) {
			 				copia.elementos.setElementAt((E) copia.elementos.elementAt(pfs).multiplicarInversa(epivot), pfs);
			 				pfs--;
			 				if (pfs < 0)
			 					break;
						}
					while(f == identidad.filas[pfsI] && identidad.filas[pfaI] == i) {

						if(identidad.columnas[pfaI] == identidad.columnas[pfsI]){
							E mult = (E) identidad.elementos.elementAt(pfsI).multiplicar(epivot);
							identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfaI).multiplicar(oper).restar(mult), pfsI);
							if (identidad.elementos.elementAt(pfsI).esCero()) {
								identidad.eliminarElem(pfsI);
								pfaI--;
								ppI--;

							}
							pfaI--;
							pfsI--;
						} else if (identidad.columnas[pfaI] < identidad.columnas[pfsI]) {
							identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfsI).multiplicarInversa(epivot), pfsI);
							pfsI--;
						} else {
							identidad.cargarElemento((E) identidad.elementos.elementAt(pfaI).multiplicar(oper), f, identidad.columnas[pfaI]);
							ppI++;

						}
						if (pfsI < 0)
							break;


					}
						if (pfsI >= 0)
						while(f == identidad.filas[pfsI]) {
						    identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pfsI).multiplicarInversa(epivot), pfsI);
							pfsI--;
							if (pfsI < 0)
								break;
						}
						while(identidad.filas[pfaI] == i) {
							identidad.cargarElemento((E) identidad.elementos.elementAt(pfaI).multiplicar(oper), f, identidad.columnas[pfaI]);
							ppI++;

						}

				} else
					while(f == copia.filas[pfs]) {
						pfs--;
						if (pfs < 0)
							break;
					}


			}

		}
		epivot = (E) copia.elementos.elementAt(0);
		int pos = 0;
		f = 1;
		while(f == 1){
			identidad.elementos.setElementAt((E) identidad.elementos.elementAt(pos).dividir(epivot), pos);
			f = identidad.filas[++pos];
		}
		return identidad;
	}
	
	/**
	 * Traspone la Matriz
	 */
	public void traspuesta() {
		// Invierte los vectores fila y columna y los reordena
		int [] aux = filas;
		filas = columnas;
		columnas = aux;
		quicksort(0, n - 1);
	}
	
	/**
	 * Devuelve el Orden de la Matriz
	 * @return Orden de la Matriz
	 */
	public int getOrden() {
		return orden;
	}

	/**
	 * Imprime la Matriz
	 */
	public void imprimirMatriz() {
		for(int i = 0; i < n; i++)
			System.out.println("" + i + "-> " + elementos.elementAt(i) + ", " + filas[i] + ", " + columnas[i]);
	}
	
	/**
	 * Devuelve el Elemento en esa posici�n
	 * @param fila del Elemento
	 * @param columna del Elemento
	 * @return Elemento de esa posici�n
	 */
	public E obtenerElemento(int fila, int columna) throws Exception {
		
		// Verifica que este dentro de los l�mites de la Matriz
		if (fila > orden || columna > orden || fila < 1 || columna < 1)
			throw new Exception("Posicion fuera de rango");
		
		E cero = elementos.elementAt(0); // Solo para crear la variable cero en memoria
		int pos = binarySearch(fila, columna); // Hace una b�squeda binaria y devuelve la posici�n del elemento. Si no lo encuentra devuelve -1
		return (pos > -1 ? elementos.elementAt(pos): (E) cero.devuelveCero()); // Si no existe en los vectores, el elemento es cero
	}

    /**
     * Method to swap to elements in an array. [Weiss]
     * @param a an array of objects.
     * @param index1 the index of the first object.
     * @param index2 the index of the second object.
     */
    private final void swapReferences(int index1, int index2 )
    {
        int tmp = filas[index1];
        filas[index1] = filas[index2];
        filas[index2] = tmp;

		tmp = columnas[index1];
        columnas[index1] = columnas[index2];
        columnas[index2] = tmp;

		E tmpe = elementos.elementAt(index1);
        elementos.setElementAt(elementos.elementAt(index2), index1);
        elementos.setElementAt(tmpe, index2);
    }

    /**
     * Return median of left, center, and right. [Weiss]
     * Order these and hide the pivot.
     */
    private int median3(int left, int right)
    {
        int center = (left + right) / 2;
        if( filas[center] < filas[left] || (filas[center] == filas[left] && columnas[center] < columnas[left]))
            swapReferences(left, center );
        if( filas[right] < filas[left] || (filas[right] == filas[left] && columnas[right] < columnas[left]))
            swapReferences( left, right );
        if( filas[right] < filas[center] || (filas[right] == filas[center] && columnas[right] < columnas[center]))
            swapReferences(center, right);
        swapReferences(center, right - 1);
        return right - 1 ;
    }

    /**
     * Internal quicksort method that makes recursive calls.
     * Uses median-of-three partitioning and a cutoff of 10. [Weiss]
     * @param left the left-most index of the subarray.
     * @param right the right-most index of the subarray.
     */
    private void quicksort(int left, int right)
    {
        if(left + CUTOFF <= right)
        {
			 int pospiv = median3(left, right);
	         int pivotfilas = filas[pospiv];
			 int pivotcolumnas = columnas[pospiv];
                // Begin partitioning
          int i = left, j = right - 1;
          for( ; ; )
            {
                while(filas[++i] < pivotfilas || (filas[i] == pivotfilas && columnas[i] < pivotcolumnas)) {}
                while(filas[--j] > pivotfilas || (filas[i] == pivotfilas && columnas[j] > pivotcolumnas)) {}
                if( i < j )
	                swapReferences(i, j);
                else
   	              break;
            }

            swapReferences(i, right - 1);   // Restore pivot
	        quicksort(left, i - 1);    // Sort small elements
	        quicksort(i + 1, right);   // Sort large elements
        }
        else  // Do an insertion sort on the subarray
	         insertionSort(left, right);
    }

    /**
     * Internal insertion sort routine for subarrays
     * that is used by quicksort. [Weiss]
     * @param a an array of Comparable items.
     * @param left the left-most index of the subarray.
     * @param right the right-most index of the subarray.
     */
    private void insertionSort(int left, int right)
    {
        for(int p = left + 1; p <= right; p++)
        {
            int tmpfilas = filas[p];
            int	tmpcolumnas = columnas[p];
            E tmpelementos = elementos.elementAt(p);
            int j;

            for(j = p; j > left && (tmpfilas < filas[j - 1] || (tmpfilas == filas[j - 1] && tmpcolumnas < columnas[ j - 1 ])) ; j--){
                filas[j] = filas[j - 1];
            	columnas[j] = columnas[j - 1];
            	elementos.setElementAt(elementos.elementAt(j - 1), j);
            }
            filas[j] = tmpfilas;
            elementos.setElementAt(tmpelementos, j);
            columnas[j] = tmpcolumnas;
        }
    }
    
    /**
     * M�todo interno para duplicar los vectores Fila, Columna y Elementos
     */
	private void duplicarVectores() {

		int [] orig_filas = filas;
		int [] orig_columnas = columnas;

		int longitud = filas.length;
		int nuevoTamano = filas.length + CAPACIDAD_POR_DEFECTO;

		filas = new int [nuevoTamano];
		columnas = new int [nuevoTamano];

		for (int i = 0; i < longitud; i++) {
			filas[i] = orig_filas[i];
			columnas[i] = orig_columnas[i];
		}
	}
	
    /**
     * M�todo interno para hacer un hueco en los vectores
     * Se utiliza para insertar elementos ordenadamente
     */
	private void huecoVectores(int posicion) {
		
		// Mueve todos los elementos una posici�n para crear un hueco
		for (int i = n; i > posicion; i--) {
			filas[i] = filas[i - 1];
			columnas[i] = columnas[i - 1];
		}

	}
	
	/**
     * M�todo interno para eliminar un elemento
     * @param Posici�n del Elemento en los vectores
     */
	private void eliminarElem(int posicion) {
		
		elementos.removeElementAt(posicion);
		
		// Para eliminar, mueve todos los elementos una posici�n sobreescribiendo el que deseamos borrar
		for (int i = posicion; i < n; i++) {
			filas[i] = filas[i + 1];
			columnas[i] = columnas[i + 1];
		}
		n--;
	}
	
	/**
     * M�todo interno binarySearch
     * Busca en los vectores columna y fila
     * @param x Limite izquierdo
     * @param y Limiete Derecho
     * @return Elemento encontrado
     */
	private int binarySearch(int x, int y) {

		int low = 0;
		int high = n - 1;
		int mid;

		while(low <= high)
		{
			mid = (low + high) / 2;

			if(filas[mid] < x)
				low = mid + 1;
			else if(filas[mid] > x)
				high = mid - 1;
			else { // Si encontr� la fila, busca la columna
				if(columnas[mid] < y)
					low = mid + 1;
				else if(columnas[mid] > y)
					high = mid - 1;
				else
					return mid;
			}
		}

		return -1;
	}

	/**
     * Verifica que en la fila exista un elemento distinto de cero
     * @param x Elemento
     * @throws No existe fila
     */
	private void existeFila(int x) throws Exception {
		int r = binarySearch(x);
		if (r == -1) // Si devolvi� -1, no existe fila
			throw new Exception ("No existe fila");
	}

	/**
     * M�todo interno binarySearch
     * Busca solo en las filas
     * @param x Elemento
     * @return Elemento encontrado
     */
	private int binarySearch(int x) {

		int low = 0;
		int high = n - 1;
		int mid;

		while(low <= high)
		{
			mid = (low + high) / 2;

			if(filas[mid] < x)
				low = mid + 1;
			else if(filas[mid] > x)
				high = mid - 1;
			else
				return mid;
		}

		return -1;
	}
	
	/**
     * M�todo interno para hacer una copia de la Matriz
     * @return Copia de la Matriz
     */
	private MatrizDispersa copiar() {
		// Crea una nueva matriz dispersa y copia cada uno de los vectores
		MatrizDispersa<E> mat = new MatrizDispersa<E>(orden);
		for(int i = 0; i < n; i++) {
			mat.elementos.insertElementAt(elementos.elementAt(i), i);
			mat.filas[i] = filas[i];
			mat.columnas[i] = columnas[i];
		}
		mat.n = n;
		return mat;
	}
	
	private int orden; // Orden de la matriz
	private int CAPACIDAD_POR_DEFECTO = 1000; // Capacidad por defecto de Vectores
	private Vector<E> elementos; // Vector de elementos
	private int []filas; // Vector de filas
	private int []columnas; // Vector de columnas
	private int n; // Cantidad de Elementos distintos de cero
    private static final int CUTOFF = 3; // L�mite del quicksort
}
